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In system analysis -work there arises from time to time a need for sequences 
of random numbers, to simulate, e.g., the effects of noise and errors of 
a random nature on system performance; and, of course, such sequences are 
essential to the use of Monte Carlo techniques. This memorandum will discuss 
some techniques for generating sequences of numbers which are suitable for 
these puproses, and will describe some experiments that have been conducted 
to study the nature of such sequences. Although the following discussion will 
be limited to application to Autone tics' RECOMP digital computer, the basic 
ideas are certainly applicable to digital computation in general. 

A conventional digital computer such as RECOMP is, of course, incapable of 
performing a truly random process; all of its operations are deterministic. 
However, there exists a rather convenient method of generating a sequence of 
numbers which, from the standpoint of the user are (through his ignorance, 
if you will) unpredictable and in this sense are pseudo-random. By a pseudo - 
random sequence we mean a previously determined sequence which is used to 
simulate a random sequence; however, in the discussion that follows the pre- 
fix "pseudo" will be omitted when we refer to such sequences. 

The method for generating this sequence is as follows: Let N be an integer 
greater than one and let x Q be a fraction 

c x Cl 



Define 

x. + _ * fractional part of 

Then the sequence /x^ j is uniformly distributed between zero and one. 



Nx k 



That is, the probability density function 
p(x) 



)1 < x < 1 



to 



otherwise 



To generate this random sequence on RECOMP we store the fraction x at 
a binary scale of zero and the integer N at a binary scale of 39. Then 
if we multiply N by x we have the fractional part of the product, which 
is the new x, in the R register. The sequence of commands is 

CLA x 

MPY N 

XAR (exchange A & R registers) 

STO x 
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The simplicity of this method is apparent. In practice it is recommended 
that X • 2~3° and that an odd power of 3 or 5 be selected for N. Different 
choices will of course provide different sequences. It is best to select 
the largest odd power of 3 or $ that can be contained in 39 bits. For further 
discussion on this point, as well as the mathematical nature of the pseudo- 
random sequence, see references (b) and (c). A convenient method of obtaining 
different sequences is to employ two generators with different odd powers of, 
say, 5 for N. A member of the first sequence is selected, with an element of 
chance, by throwing a sense switch, to be used as the starting number for the 
second sequence which provides the random numbers for the problem. A brief 
table of odd powers of 3 and 5 in command format is given in Appendix fi. 

In a practical application we are concerned with two characteristics of the 
random sequence. One of these is the distribution of the random numbers. As 
asserted above the members of the random sequence are uniformly distributed 
between zero and one, and a proof of this fact may be found in reference (d). 
Figure 1 shows the actual distribution of a sample of 102U numbers generated 
in this manner. 

A second characteristic of interest is the sequencing of the numbers, or more 
precisely, if the sequence is thought of as a random time series, the power 
spectral density of the series. In this regard experimental evidence indicates 
that the spectral density is "white" or uniform over all frequencies. (Of 
course, as in all digital computer work, the spectral content is band limited 
in a real time sense, by the sampling frequency, or the rate at which the 
numbers are generated.) Another way of considering this characteristic which 
does not depend on any concept of time is to state that the members of the 
sequence are statistically independent of one another. To demonstrate this 
fact an "autocorrelation" function was computed. 

N 

H ( ) - i K ■ *k V -012... 
x N K ■ 1 

Figure 2 shows an actual example of this computation on a sample of 102U 
numbers. (The numbers were shifted first by subtracting one-half to lie 
between plus and minus one-half.) From this figure it is apparent that for 
one or more shifts the numbers are uncorrelated. 

An important consequence of the independence of the members of the sequence 
is that it permits the use of a single random number generator to provide 
numbers for several applications in the same problem. 
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Figure 1. Sample Distribution of 102U Random Nutabers 
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Thus, we have available a simple method of generating a random sequence whose 
members are uniformly distributed between zero and one and are independent of 
one another or "white". Other rectangular distributions may be easily obtained 
by multiplying by a scale factor and adding a bias. For example, suppose it 
is desired to select at random an integer between 1 anr) $2, We simply generate 
a random number, multiply by 52, add 1, and take the integral part of the answer 
as the desired number. 

However, distributions other than rectangular are often required. For example, 
failure rate is often characterized by an exponential distribution, system 
errors are frequently considered to have gaussian distributions; other random 
events may have a Poisson distribution. How may other distributions be ob- 
tained from a rectangular distribution? 

If x is uniformly distributed between zero and one and if the probability 
density function p (t) has the properties 



p (t) > and I p (t) dt = 1 



then the random variable z defined by the relation 

z 



x * I p(t) dt (1) 



has the probability distribution function 

z 
P (z") =1 p (t) dt 



and hence the probability density function p (z). For Probability >z < zt 



Probability f x < f p (t ) dt 



; 



z o 



p (t) dt 



since x is uniformly distributed between zero and one and hence Probability 
ix •£- ac s a 



For example, suppose an exponential distribution 


is required, i.e., 


C -a. 




,----. \ae 
p(z) » <) 




& z 

z <. 


V 
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From equation (1) it is therefore required that the uniform random variable 

ae~ at dt " . 

£ z 

z <C 

\ 

For z £ 

t -az 
x ■ 1 - e 

or z ■ 1 log (l-x) 

-a 

or since x and (1 - x) have the same distribution, -we let 

z ■ 1 log x 
-a 

and the random variable z will have the probability density p (z) as 
desired. In other words, we generate a random number, compute its natural 
logarithm, and divide by minus a • The result will have the probability 
density function p (z). 

As a second example suppose it is desired to select a point at random from 
the unit circle under the assumption that the points are uniformly distributed 
over the circle. One possible solution would be to generate two random numbers, 
say x and y, scaled and biased to lie between plus and minus one. rejecting 
the pair if the sum of their squares exceeded one. The pairs (x, y) that were 
accepted would, of course, have a uniform distribution over the unit circle* 

An alternate method would use polar coordinates r and 9. From symmetry it 
is clear that any angle is equally likely so the probability density of & 
is given by 



'*»> 



1 0<8 £ 2TT 

2TT 

otherwise 



Therefore, to generate the random variable Qwe simply generate a random 
number and multiply it by 2jT* 
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To find the distribution of r we note that the probability of a point lying 
in an incremental area is 



dx dy 


B rdrd^ 


IT 


TT 




* (2r dr) d9 




2TT 




* p (r)dr p (0) de 



Since we have already determined p^( G) it follows that 

2r < r -1. 1 




otherwise 

Prom equation (l) we require that the uniformly distributed number 

r 
x - | 2t dt = r 2 
^0 

or r * x 

To summarize the procedure then we generate two random numbers. One of these 
is multiplied by 2 TT and designated Q . The square root of the other is 
extracted and the result designated r. The resulting points (r, B) are 
uniformly distributed over the unit circle. The size of the circle can be 
easily scaled by multiplying r by the radius of the desired circle. 

In principle any desired distribution maybe obtained by the method discussed 
in the preceding paragraphs. However, there may be computational difficulties 
for some density functions. Unfortunately, this is the case for the most 
important gaussian distribution, 

1 -2/2 
p(z) - e 

According to equation (1) it would be necessary to invert the equation 

C -t 2 / 2 

x * I e dt 
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If the function on the right is designated 
to find 



(z) we see that it is required 



-1 

(!) (x) 



and of course this equation may be solved by numerical methods? but, it would 
appear that considerable computation would be required, and for this reason it 
has not been attempted. As an alternate approach we again resort to polar 
coordinates. If u and v are independent random numbers with gaussian 
distributions then the random variable 

\T"2 ♦ 2 1 

r - >|u v 

has the so-called Rayleigh distribution 

fre~ r A r >i 

p (r) - i 

r -c 



and the random variable 



-1 



u 



is uniformly distributed between zero and 2 IT • To obtain a random number 
with the Rayleigh distribution, as before, we let the uniformly distributed 
number 

i ( 2 



te 



dt 



- e-% 



or 



2ioe y=) 

or since x and (1-x) have the same distribution we simply let 



i 
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I 



1 



Similarly we let 

a - s 

where of course x is another uniformly distributed number. Then we simply 
make the inverse transformation 



u = r cos Q 

v ■ r sin Q 

and u and v will be independent random numbers with gaussian distributions. 
It will be noted that this method requires the computation of a logarithm and 
sine and cosine for which subroutines are normally available. 

The above method provides a pair of independent gaussian random numbers at the 
expense of computing a logarithm and sine-cosine. Since the latter computations 
are somewhat time consuming, and usually time is at a premium in analyses of 
a probabilistic nature, this method has limited practical application. The 
following approach provides a technique for generating numbers whose distribution 
is approximately gaussian with relatively little computation required. 

From the central limit theorem of statistics it is known that the distribution 
°^ *^ e av ®rage of N samples from any distribution approaches the gaussian 
distribution as N becomes large. If the original distribution is rectangular, 
the convergence is remarkably fast, in fact, N * 3 or k gives a distribution 
which is a very good approximation to the gaussian. 



Let 



1 04-x^.l 



Pi (x-t ) - -j (2) 

otherwise 



for i - 1, 2, . . . , N 
Then the random variable 



N 



x 



1 
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has the probability density function p n (z) which may be obtained by con- 
volving p 1 with p^ or 



p (z) - 
n 

J 



p (z-x) p (x) dx N - 2, 3, . >\. 
90 



For convenience we let 

tt ' z -^n 

where ,// n " J tP n (*)dt 






(t "An )2 P n (t)dt 



so that il has a zero mean and unit variance. 

For N ■ 2 the distribution is triangular shaped and of limited interest. For 
N ■ 3 evaluation of the convolution integral gives 

3 - u 2 < |u|< 3 

8 "" 'T 

P (u) « \ (3- [uf ) 2 14 N(£ 3 (3) 

16 

3^ a /u'/ 

where u * 2 (x + x + x ) - 3 (it) 

with x., i ■ 1> 2, .3* distributed according to (2). A plot of p^ (u) is 
given in Figure 3. For comparison, a plot of the gaussian distribution 

is given on the same figure. From this figure it is seen that the distribution 
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Figure 3» Comparison of Distribution p^ (u) with Gaussian, 
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of u is a very good approximation to the gaussian. In fact, it would be 
difficult to distinguish between samples from the two distributions unless 
an extremely large sample were taken. Figure k shows the actual distri- 
bution of ^000 numbers generated according to equation (li). In this 
figure the area of the rectangles equals the fraction of the sample falling 
in the corresponding interval. The gaussian is also plotted on this figure 
for comparison. 

This technique lends itself readily to digital computation. A subroutine for 
generating "approximate gaussian" random numbers in this fashion is given in 
Appendix A. As mentioned above, these numbers have the distribution given 
by (3) with zero mean and unit variance. To simulate the gaussian random 



variable z with mean JL( and variance 






let 



z - CTU + M„ 
z z 

It should be noted that members of the resulting gaussian random sequence 
are also independent of one another or "white". 

An even better approximation may be obtained by combining four of the uni- 
formly distributed numbers. For N equal U, evaluation of the convolution 
integral yields 



p u 



(u) 



^ h ff- 2 p 1 u +/u/ 



</ (2 l7- /W ) 3 
5u 



OS juj< *yT 



where u ■ 






h 

V""l 

T"l 1 



- 2 



with x^, i ■ 1, 2, 3, U, distributed according to (2). As before the distri- 
bution has been scaled to have a zero mean and unit variance. Figure £ shows 
a plot of pi (u) with the gaussian for comparison. 

Appendix C provides a tabulation of the density functions p~ (u) and p, (u) 
together with their respective cumulative distributions. 
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Another technique that is useful in providing a distribution that, although 
not gaussian, favors small numbers over large and thus may be adequate for 
some purposes, is to simply multiDly two of the uniformly distributed numbers. 
If x,, Xp are distributed according to 





U 


-1 ^,X. -_ 


P(x ± ) « 


i 


i 




c° 


otherwise 


i - 1, 2 





then the random variable z * x v 

1 2 



has the distribution 



p(z) 



h log J^ -1 ... z . . 1 

izj 



\ otherwise 

Also the random variable y = x^ j x-, | has the distribution 

P (y) ■ 



T^T-h -i^y-i 



/ otherwise 

\ 

Random variables such as these have the advantage of being computed rather 
easily. 

As a final subject, we will briefly consider the case where a random sequence 
is required with a spectral density other than white. For example, it might 
be required to constrain the sequence so that it does not change value too 
rapidly or in other words, suppress the high frequency components in the sequence. 
To achieve this end it is necessary to run the "white" sequence through a low- 
pass filter. 

If the input x (t) to a filter with transfer function H (jw) is white, i.e., 
has power spectral density, 

S (jw) » S (0) 
x x 

then the power spectral density of the output y (t) is 

S (Jw) -Jh (jw) I 2 S (0) 
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or in other words, the output spectral content is determined by the filter 
characteristic • If the input to a linear filter has a gaussian distribution 
then the output will also be gaussian. It is then only necessary to determine 
the mean and variance of the output in order to completely characterize the 
output distribution. Well known techniques are available for the design of 
digital filters and thus it is possible to generate gaussian random numbers 
with a desired spectral content. 

To demonstrate this technique a simple first order lag filter was programmed 
and fed with a white gaussian input. The distribution and autocorrelation 
function of the output were then computed. Such a filter has the character- 
istic 

H(jw) = a 

a + jw 



and impulse response 



h(t) = ae~ 



where a is the so-called corner frequency. The output spectral density 
is therefore 

s y (jw) s a2 S x (0) (7) 

2 2 

a^+ w* 1 

The autocorrelation function of the output(the inverse transform of S (jw) 

y 

has the form _ , , n N 

, s 2 -an*/ (8) 

V £> - v 7 e ' 

where {Y is the variance of the output. 

To derive the difference equation defining the digital filter, we use the 
fact that the output is the convolution of the input with the filter impulse 
response, 

-t 



y(t) * ( h(t - t«) x(t') dt» 

Jo 



&. 



\ t ae" a(t "" t,} x(t') dt' 
JO 
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To evaluate this integral numerically we take an integration step A t s T 
and assume x (t) is constant over this interval, thus 

x(t) - x n (n-l)T -'. t < nT n * 1, 2, . .. 

The x will be the input white gaussian random sequence, with zero mean and 
unit variance. 

Now 

(n+l)T 

, r x . -a(n+l)T f at 1 , . 

y n4 .i 3 y(fe + l)T) = e \ ae x(t»)dt' 

J 



nT f n+1 )t 

-at ^-anTf ^at' / x ^ ja ., . -a(n+l)Tf V 'at' 



-ax. -arm ax,' /.,\,., ^ -avn+in ax- 
= e e i ae x (t')dt' + e ' { ae x^^dt 1 



JnT 



or 



-aT . /^ -aT< 



^n+1 - e" al y n «■ (1 - e ~ ai ) x 



n+1 

Equation (9) is the difference equation that defines the filter. By averaging 
both sides of (9) we see that the output mean 

or /-'( « /-I ■ since the input has zero mean. The variance of the output 

y x 

may be obtained by squaring both sides and averaging, and noting that the cross 
terra on the right involving x y averages to zero, since future input x + -> 

is uncorrelated with present output y . 



Thus 



CT 2 = e - 2aT 2 v ♦ (1 - e" aT ) 2 cT * 



o t „-aT 2 

or _2 _ 1 - e __ 

u y 1 + e -aT ' x 
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Prom equations (8) and (10) the autocorrelation function of the output is 



t -aT 
1 - e 


e -raT 


y 1 + e" aT 




for H, - 0, 1, 2, . . . 




As an example of this computation we let 




aT - 0.1* 




This gives 




2 




j- * .197 




y 





Figure 6 shows the distribution of the output (normalized to have unit variance) 
and Figure 7 shows the output autocorrelation function. One could apply Fourier 
techniques to determine the output spectral content. However, from the exponen- 
tial nature of the output autocorrelation function, it is clear that the output 
has the desired spectral density. 
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Figure 6, Output Distribution, Filtered Gaussian Random Numbers. 
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Figure 7. Output Autocorrelation Function, Filtered Gaussian Random Numbers. 



o 



t-3 



RECOMP TECHNICAL BULLETIN NO. 22 PAGE TWENTY-ONE 



Appendix A "Approximate - Gaussian" Random Number 

Generator Subroutine 



Enter Subroutine with ♦ TRA 0006.0 

Exit to next location with gaussian random number in A and R registers 



0006.0 
0010.0 



0020.0 



+ 


SAX 7760.0 


+ 


CTL 


0010.0 


+ 


CTV 0020.0 


+ 


TRA 


7760.0 


+ 


ADD 7773.0 


+ 


STA 


7773.1 


+ 


CLA 7777.0 


+ 


MPY 


7776.0 


+ 


XAR 0000.0 


+ 


STO 


7777.0 


+ 


ARS 0001.0 


+ 


ADD 


777h.O 


+ 


STO 777U.O 


•f ' 


CLA 


7777.0 


+ 


MFY 7776.0 


+ 


XAR 


0000.0 


+ 


.STO 7777.0 


+ 


ARS 


0001.0 


+ 


ADD 777li.O 


+ 


STO 


777^.0 


+ 


CLA 7777.0 


+ 


MPY 


7776.0 


•f 


CU 777£.0 


+ 


XAR 


0000.0 


+ 


STO 0027.0 


•f 


ARS 


0001.0 


•f 


FAD 777Ii. 


+ 


TRA 


0000.1 



-3 at Binary scale of 2 
+2 at Binary scale of 39 
N at Binary scale of 39 
x Q ■ +1 at Binary scale of 39 



Notes J 

(l) It is recommended that N be an odd power of 3 or J> # 



(2) Main memory addresses are underlined. 
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Appendix E 



Odd Powers of 3 and 5 



RECOMP users may find the following list of odd powers of 3 and 5> in 
Command formats useful in programming random number generators: 



3 23 


■ «f 


1275321 


+ 


6720U51 


3 21 


m + 


0115731 


+ 


U2U5311 


3 19 


m + 


0010520 


- 


651*7551 


3 17 


■ + 


0000751 


- 


2U13U11 


3 1S 


■ + 


0000061 


♦ 


2714*651 


5 15 


« + 


03h3271 


+ 


5223061 


P 


■ + 


0011060 


*■ 


2317121 


JX 


a + 


0000270 


+ 


1035561 


S 9 


■ + 


oooooca 


+ 


5632621 
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Appendix C Probability Density Function and Cumulative 

Distribution of p^(u) and Pk(u) 









u^ 








,..L 


u 


P 3 (u 


l) 


\ p^Mdx 

o 




v h M 






.0 


+.37500 


+ 


+.00000 + 





+ .6,6661 + 





+.00000 + 


.1 


+.37375 


+ 


+ .371*58 - 


1 


+.6631*3 + 





+ .381x27 - 1 


.2 


+.37000 


+ 


+ .71x667 - 


1 


+.651*10 + 





+ .761*89 - 1 


.3 


+.36375 


+ 


+.11137 + 





+.63926 + 





+.11385 + 


.1* 


+.35500 


+ 


+ .11*733 + 





+ .6191*9 + 





+.15021 + 


.5 


+.31*375 


+ 


+.18229 + 





+.59536 + 





+.18530 + 


.6 


+.33000 


+ 


+.21600 + 





+ .5671*5 + 





+.21888 + 


.7 


+.31375 


+ Q 


+ .21*821 + 





+.536311 + 





+.25076 + 


.8 


+.29500 


+ 


+.27867 + 





+.50260 + 





+.23076 + 


.9 


+.27375 


+ 


+.30712 + 





+.1*6681 + 





+ .308 76 + '.) 


1.0 


+.25000 


+ 


+.33333 + 





+ .1*2956 + 





+.331*61* + 


1.1 


+.22563 


+ 


+.35710 + 





+ .3911*1 + 





+ .35831* + 


1.2 


+.20250 


+ 


+.37850 + 





+.35291* + 





+.37983 + 


1.3 


+.18063 


+ 


+.39765 + 





+.311*71* + 





+.39910 + 


1.1* 


+.16000 


+ 


+ .I1II1 67 + 





+.27737 + 





+ .1*1619 + 


l.S 


+ .11*063 


+ 


+ .1*2969 + 





+ .2lili*3 + 





+.1*3116 + 


1.6 


+.12250 


+ 


+ .lil*283 + 





+ .2071*7 + 





+ .1*1*1*10 + 


1.7 


+.10563 


+ 


+ .ii51*23 + 





+.17609 + 





+.15^16 + 


1.8 


+.90000 


- 1 


+ . I16I4.OO + 





+. lb 781 + 





+ .1*61*50 + 


1.9 


+.75625 


- 1 


+ .1*7227 + 





+.12273 + 





+ .1*7229 + 


2.0 


+.62500 


- 1 


+.1*7917 + 





+.10067 + 





+ .1*7873 + 


2.1 


+.50625 


- 1 


+.l*81i8l + 





+ .811*15 - 


1 


+ .1*8397 + 


2.2 


+.1*0000 


- 1 


+ .1*8933 + 





+ .61*791 - 


1 


+ .1*8818 + 


2.3 


+.30625 


- 1 


+ .1*9285 + 





+.50599 - 


1 


+ .1*9150 + 


2.1* 


+.22500 


- 1 


+ .1*9550 + 





+ .3861*7 - 


1 


+ .1*91*06 + 


2.5 


+.15625 


- 1 


+ .1*971*0 + 





+ .2871*3 - 


1 


+ .1*9600 + 


2.6 


+.10000 


- 1 


+ .1*9867 + 





+.20695 - 


1 


+ .1*971*2 + 


2.7 


+.56250 


- 2 


+ .1*991x1* + 





+.11*309 - 


1 


+ .1*981*2 + 


2.8 


+.25000 


- 2 


+ .1*9983 + 





+.9391*1* - 


2 


+.L9910 + 


2.9 


+.62500 


- 3 


+ .1*9998 + 





+.57576 - 


2 


+ .1*9953 + 


3.0 


+.81*703 


-21 


+.50000 + 





+.32063 - 


2 


+ .1*9979 + 


3.1 


+.00000 


+ 


+.10000 + 


1 


+ .151*82 - 


2 


+.1*9992 + 


3.2 


+.00000 


+ 


+.10000 + 


1 


+.59085 - 


3 


+ .1*9998 + 


3.3 


+.00000 


+ 


+.10000 + 


1 


+.llil7l* - 


3 


+.50000 + 


3.1 


+.00000 


+ 


+.10000 + 


1 


+.81*1*81* - 


5 


+.50000 + 


3.5 


+.00000 


+ 


+.10000 + 


1 


+.00000 + 





+.50000 + 


3.6 


+.00000 


+ 


+.10000 + 


1 


+.00000 + 





+.50000 + 
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